In-class Exercise 5

1. Install packages

packages <- c('maptools', 'sf', 'raster', 'spatstat', 'tmap', 'tidyverse','plotly', 'ggthemes')
for(p in packages){
  if(!require(p, character.only=T)){
    install.packages(p)
  }
library(p,character.only = T)
}

• Note:
  • Maptool is for creating spatical objects
  • sf is for importing geospatial data
  • raster is for creating raster object
  • spatstat is for performing SPPA
  • tmap is for plotting thematic maps
  • tidyverse is for data handling

2. Import data

CHAS <- read_rds("data/rds/CHAS.rds")
childcare <- read_rds("data/rds/childcare.rds")

mpsz_sf <- st_read(dsn = "data/shapefile", layer = "MP14_SUBZONE_WEB_PL")

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\yiling-yu\IS415_Blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

glimpse(CHAS)

Rows: 1,192
Columns: 17
$ NAME                    <chr> "FONG CLINIC", "FRASER MEDICAL CENTR~
$ DESCRIPTION             <chr> "CDMP,CHAS", "CDMP,CHAS", "CDMP,CHAS~
$ ADDRESSBLOCKHOUSENUMBER <chr> "162", "78A", "1", "982", "163", "30~
$ ADDRESSFLOORNUMBER      <chr> "01", "01", "01", "01", "01", "01", ~
$ ADDRESSPOSTALCODE       <chr> "140162", "101078", "150001", "53098~
$ ADDRESSSTREETNAME       <chr> "MEI LING ST", "TELOK BLANGAH ST 32"~
$ ADDRESSUNITNUMBER       <chr> "361", "07", "4524", "03", "426", "8~
$ X_COORDINATE            <chr> "24609.88367388", "25288.1926698", "~
$ Y_COORDINATE            <chr> "30516.49536544", "28423.68094962", ~
$ HCI_CODE                <chr> "9403524", "16C0164", "9400550", "94~
$ LICENCE_TYPE            <chr> "MC", "MC", "MC", "MC", "MC", "MC", ~
$ HCI_TEL                 <chr> "64745993", "62733603", "62726628", ~
$ ADDR_TYPE               <chr> "A", "A", "A", "A", "A", "A", "A", "~
$ Lat                     <chr> "1.2922546717017", "1.2733280593593"~
$ Lng                     <chr> "103.802856952292", "103.80895202390~
$ ICON_NAME               <chr> "OM2_MOHicon_CHASclinic_20px.png", "~
$ ADDRESSBUILDINGNAME     <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, ~

glimpse(childcare)

Rows: 1,545
Columns: 7
$ NAME              <chr> "AVERBEL CHILD DEVELOPMENT CENTRE PTE LTD"~
$ DESCRIPTION       <chr> "Child Care Services", "Child Care Service~
$ ADDRESSPOSTALCODE <chr> "760742", "159053", "556912", "569139", "4~
$ ADDRESSSTREETNAME <chr> "742, YISHUN AVENUE 5, #01 - 470, SINGAPOR~
$ Lat               <chr> "1.42972006815634", "1.28668026511187", "1~
$ Lng               <chr> "103.833109448847", "103.813766346615", "1~
$ ICON_NAME         <chr> "onemap-fc-childcare.png", "onemap-fc-chil~

mpsz_sf <- st_set_crs(mpsz_sf, 3414)
st_crs(mpsz_sf)

Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

• WGS84 for childcare dataset is assumed. why?

3. Coverting the aspatial data frame into sf objects

CHAS_sf <- st_as_sf(CHAS,
                    coords = c("X_COORDINATE",
                               "Y_COORDINATE"),
                    crs=3414)

• Note: st_as_sf() accept coordinate in character data type

childcare_sf <- st_as_sf(childcare,
                    coords = c("Lng",
                               "Lat"),
                    crs=4326) %>%
  st_transform(crs = 3414)

#plot the graph to have a look
tmap_mode("view")
tm_shape(childcare_sf) +
  tm_dots(alpha=0.4,
          col= "blue",
          size = 0.05) +
tm_shape(CHAS_sf) +
  tm_dots(alpha=0.4,
          col= "red",
          size = 0.05)

4. Convert sf to ppp objects

a. Convert sf data frame to Spatial class

childcare <- as_Spatial(childcare_sf)
CHAS <- as_Spatial(CHAS_sf)
mpsz <- as_Spatial(mpsz_sf)

b. Convert Spatial class into generic sp object

childcare_sp <- as(childcare, "SpatialPoints")
CHAS_sp <- as(CHAS, "SpatialPoints")
mpsz_sp <- as(mpsz, "SpatialPolygons")

• as.SpatialPoint() of as.SpatialPolygon() of maptools package

c. Covert generic sp object into ppp object

childcare_ppp <- as(childcare_sp,"ppp")
CHAS_ppp <- as(CHAS_sp,"ppp")

• using as.ppp() of maptools package.

d. Remove duplicate points using jitter

childcare_ppp_jit <- rjitter(childcare_ppp,
                             retry = TRUE,
                             nsim = 1,
                             drop = TRUE)
any(duplicated(childcare_ppp_jit))

[1] FALSE

CHAS_ppp_jit <- rjitter(CHAS_ppp,
                             retry = TRUE,
                             nsim = 1,
                             drop = TRUE)
any(duplicated(CHAS_ppp_jit))

[1] FALSE

5. Extracting Punggol Planning Area

pg <- mpsz[mpsz@data$PLN_AREA_N=="PUNGGOL",]

6. Creating owin object

a. Convert SpatialPolygonsDataFrame into SpatialPolygons object

pg_sp <- as(pg, "SpatialPolygons")

b. Convert SpatialPolygons into owin object

pg_owin <- as(pg_sp, "owin")

c. Extract spatial points window owin

childcare_pg <- childcare_ppp_jit[pg_owin]
CHAS_pg <- CHAS_ppp_jit[pg_owin]

plot(childcare_ppp_jit)

plot(childcare_pg)

7. Analysing Spatial Point Process Using L-Function

Monte Carlo test with L-fucntion

L_childcare <- envelope(childcare_pg,
                        Lest,
                        nsim = 99,
                        rank = 1,
                        gloval = TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.

Done.

L_CHAS <- envelope(CHAS_pg,
                   Lest,
                   nsim = 99,
                   rank = 1,
                   gloval = TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.

Done.

plotting interactive L-function

L_childcare

title <- "Pairwise Distance: L function"

Lcsr_df <- as.data.frame(L_childcare)

colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
  # plot observed value
  geom_line(colour=c("#4d4d4d"))+
  geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
  # plot simulation envelopes
  geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
  xlab("Distance r (m)") +
  ylab("L(r)-r") +
  geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1])  +
  geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
  geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
  theme_tufte()+
  ggtitle(title)

text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"

# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){ 
  if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text3, traces = 4) %>%
      rangeslider() 
  }else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      rangeslider() 
  }else {
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider() 
  }
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
  if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      rangeslider() 
  } else{
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider()
  }
} else{
  ggplotly(csr_plot, dynamicTicks=T) %>%
    style(text = text1, traces = 4) %>%
    style(text = text2, traces = 5) %>%
    style(text = text3, traces = 6) %>%
    rangeslider()
  }

L_CHAS

title <- "Pairwise Distance: L function"

Lcsr_df <- as.data.frame(L_CHAS)

colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
  # plot observed value
  geom_line(colour=c("#4d4d4d"))+
  geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
  # plot simulation envelopes
  geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
  xlab("Distance r (m)") +
  ylab("L(r)-r") +
  geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1])  +
  geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
  geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
  theme_tufte()+
  ggtitle(title)

text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"

# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){ 
  if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text3, traces = 4) %>%
      rangeslider() 
  }else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      rangeslider() 
  }else {
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider() 
  }
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
  if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      rangeslider() 
  } else{
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider()
  }
} else{
  ggplotly(csr_plot, dynamicTicks=T) %>%
    style(text = text1, traces = 4) %>%
    style(text = text2, traces = 5) %>%
    style(text = text3, traces = 6) %>%
    rangeslider()
}